TECHNICAL 

library 


AD  Rtto  S66 


TECHNICAL  REPORT  ARBRL-TR-02431 


FLAT  FLAME  OLYMPICS:  TEST  PROBLEM  A 


Terence  P.  Coffee 


October  1982 


US  ARMY  ARMAMENT  RESEARCH  AND  DEVELOPMENT  COMMAND 
BALLISTIC  RESEARCH  LABORATORY 

ABERDEEN  PROVING  GROUND,  MARYLAND 


Approved  for  public  release;  distribution  unlimited. 


Destroy  this  report  when  it  is  no  longer  needed. 
Do  not  return  it  to  the  originator. 


Secondary  distribution  of  this  report  is  prohibited. 


Additional  copies  of  this  report  may  be  obtained 
from  the  National  Technical  Information  Service, 

U.  S.  Department  of  Commerce,  Springfield,  Virginia 

22161. 


The  findings  in  this  report  are  not  to  be  construed  as 
an  official  Department  of  the  Army  position,  unless 
so  designated  by  other  authorized  documents. 


use  of  trade  names  or  manufacturers '  names  in  this  report 
does  not  constitute  indorsement  of  any  comteroial  product. 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER 

Technical  Report  ARBRL-TR-02431 

2.  GOVT  ACCESSION  NO. 

3.  RECIPIENT’S  CATALOG  NUMBER 

i  4.  TITLE  (and  Subtitle) 

FLAT  FLAME  OLYMPICS:  TEST  PROBLEM  A. 

S.  TYPE  OF  REPORT  ft  PERIOD  COVERED 

Final 

6.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHORfeJ 

TERENCE  P.  COFFEE 

8.  CONTRACT  OR  GRANT  NUMBERfa) 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  DRDAR-BLI 

ALpt* Ha <=»n  Proving  (TrnnnH^  MD   .  

10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  ft  WORK  UNIT  NUMBERS 

1L161102AH43 

11.  CONTROLLING  OFFICE'NAME  AND  ADDRESS 

US  Army  Armament  Research  &  Development  Command 
'  US 'Army  Ballistic  Research  Laboratory  (DRDAR-BL) 
Aberdeen  Provine  Ground.  MD  21005 

12.  REPORT  DATE 

October  1982 

13.  NUMBER  OF  PAGES 

46 

14.  MONITORING  AGENCY  N'AME  ft  ADDftESSfi/  different  from  Controffing  Office) 

IS.  SECURITY  CLASS,  (ol  thle  report) 

UNCLASSIFIED 

15a.  DECLASSIFY  CATION/ DOWN  GRADING 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (of  thle  Report) 

Approved  for  public  release;  distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  tn  Btock  20,  ft  different  from  Report) 

18.  SUPPLEMENTARY  NOTES 

19.  KEY  WORDS  (Continue  on  ravaroa  aide  if  necaeaary  and  identify  by  block  number) 

Laminar  Flame  One  Dimensional 

Premixed  Flame  Transient  Flame 

Lewis  Number  Asymptotic  Analysis 

Finite  Elements  Method  Temperature  Profiles 

20.  ABSTRACT  (XJotrtbxue  an  re  Ter  me  alxfim  ft  nee  me  maty  and  identify  by  bfock  number) 

raj 

This  report  discusses  a  test  problem  for  a  computer  program  for  numeri¬ 
cally  solving  the  equations  governing  a  laminar,  premixed,  one-dimensional 
flame.  The  problem  was  proposed  by  GAMM  (Committee  for  Numerical  Methods 
in  Fluid  Mechanics),  and  has  been  solved  for  presentation  at  a  workshop  at 
the  Technical  University,  Aachen,  Germany,  12-14  Oct.  1981. 

The  test  problem  is  an  unsteady  propagating  flame  with  one-step  chemis¬ 
try  and  Lewis  number  different  from  unity.  A  code  developed  for  steady 

DD 


FORM 
1  JAN  73 


1473 


EDITION  OF  1  NOV  65  IS  OBSOLETE 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


SECURITY  CLASSIFICATION  OF  THIS  PAGEOWifi  Pete  Entered) 


Abstract  (Cont'd): 

20.  state  problems  with  elementary  chemistry  was  modified  to  use  the 
simplified  transport  and  chemistry  of  the  test  problem  and  to  follow  the 
details  of  the  transient  solution. 


The  problem  is  solved  for  six  cases.  The  cases  differ  in  the  Lewis 


number  chosen  and  the  activation  energy  of  the  single  reaction.  The  initi 
conditions  used  are  the  steady  state  solutions  predicted  by  the  simplified 
analytic  method  of  asymptotic  analysis. 


il 


In  most  cases,  the  numerical  solutions  rapidly  converge,  and  the  stea  ly 
state  solutions  are  similar  to  the  asymptotic  solutions.  However,  in  one 
case,  with  high  activation  energy  and  Lewis  number  equal  to  two,  the  solution 
does  not  converge.  Instead,  large  oscillations  in  the  flame  speed  and  the 
profiles  occur. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE(TWi»n  Date  Entered) 


TABLE  OF  CONTENTS 


PAGE 


LIST  OF  ILLUSTRATIONS . 5 

I.  INTRODUCTION . 7 

II.  THE  TEST  PROBLEM . 8 

III.  THE  NUMERICAL  SCHEME . 11 

IV.  NUMERICAL  SOLUTIONS . 14 

V.  CONCLUSIONS . 17 

REFERENCES . 40 

DISTRIBUTION  LIST . 41 


3 


LIST  OF  ILLUSTRATIONS 


FIGURE  PAGE 

1.  CASE  1.  FLAME  SPEED  PROFILE . 18 

2.  CASE  1.  TEMPERATURE  PROFILE . 19 

3.  CASE  2.  FLAME  SPEED  PROFILE . 20 

4.  CASE  2.  FLAME  SPEED  PROFILE . 21 

5.  CASE  2.  REACTANT  PROFILE . 22 

6.  CASE  2.  TEMPERATURE  PROFILE . 23 

7.  CASE  2.  ENTHALPY  PROFILE . 24 

8.  CASE  3.  FLAME  SPEED  PROFILE . 25 

9.  CASE  3.  REACTANT  PROFILE . 26 

10.  CASE  3.  TEMPERATURE  PROFILE . 27 

11.  CASE  3.  ENTHALPY  PROFILE . 28 

12.  CASE  4.  FLAME  SPEED  PROFILE . 29 

13.  CASE  4.  TEMPERATURE  PROFILE . 30 

14.  CASE  5.  FLAME  SPEED  PROFILE . 31 

15.  CASE  5.  FLAME  SPEED  PROFILE . 32 

16.  CASE  5.  REACTANT  PROFILE . 33 

17.  CASE  5.  TEMPERATURE  PROFILE  . 34 

18.  CASE  5.  ENTHALPY  PROFILE . 35 

19.  CASE  6.  FLAME  SPEED  PROFILE . 36 

20.  CASE  6.  REACTANT  PROFILE . 37 

21.  CASE  6.  TEMPERATURE  PROFILE . 38 

22.  CASE  6.  ENTHALPY  PROFILE . 39 


5 


I.  INTRODUCTION 


We  are  involved  in  the  modeling  of  one-dimensional ,  premixed,  laminar, 
steady-state  flames.  '  Our  purpose  is  to  determine  validated  sets  of 
elementary  chemical  reactions  for  use  in  predictive  combustion  models.  This 
involves  developing  numerical  methods  to  solve  premixed  flame  problems  with 
complex  chemistry  and  transport  properties . ^ ' 4 

This  report  discusses  a  test  problem  for  the  numerical  methods.  The 
problem  was  proposed  by  GAMM  (Committee  for  Numerical  Methods  in  Fluid 
Mechanics) ,  and  has  been  solved  for  presentation  at  a  workshop  at  the 
Technical  University,  Aachen,  Germany,  12-14  Oct  81. 

The  purpose  of  the  workshop  is  to  bring  together  the  small  group  of 
scientists  actively  working  in  the  field  of  numerical  methods  in  flame 
propagation.  The  workshop  will  focus  on  numerical  methods  that  have  been 
developed  to  solve  premixed  flame  problems.  In  particular,  the  meeting  will 
seek  to  establish  the  difficulties  that  result  from  non-equal  dif fusivities  of 
heat  and  matter  (Lewis  number  not  equal  to  one)  and  the  consequences  that 
these  have  upon  the  accuracy  of  the  solution. 

At  the  workshop,  the  solution  of  two  test  problems  by  different  numerical 
schemes  will  be  compared.  Here  we  discuss  the  first  test  problem;  an  unsteady 
propagating  flame  with  one-step  chemistry  and  Lewis  number  different  from 
unity.  The  second  test  problem,  a  steady,  stoichiometric  hydrogen-air  flame 
with  complex  chemistry,  is  reported  on  elsewhere.5  The  results  of  the 
comparisons  will  be  published  in  the  proceedings  of  the  workshop.5 


2 

J .M .  Hevmevl  and  T.P.  Coffee,  "The  Detailed  Modeling  of  Premixed,  Laminar, 
Steady-State  Flames,  I.  Ozone,"  Combustion  and  Flame .  Vol.  39,  pp.  301-315, 
1980. 


2 

T.P.  Coffee  and  J.M.  Heimerl,  "Transport  Algorithms  for  Premixed,  Laminar, 
Steady  Flames Combustion  and  Flame,  Vol.  43,  pp.  273-289,  1981. 

7 

T.P.  Coffee  and  J.M  .  Heimerl,  "A  Method  for  Computing  the  Flame  Speed  for  a 
Laminar,  Premixed,  One  Dimensional  Flame,"  BEL  Technical  Report,  ARBRL-TR- 
02212,  January  1980,  (AD-A082803) . 

^ T.P .  Coffee,  "A  Computer  Code  for  the  Solution  of  the  Equations  Governing  a 
Laminar,  Premixed,  One  Dimensional  Flame,"  BRL  Memorandum  Report,  A RBRL -MR- 
031  65,  April  1982.  (AD  A114041) 

^ J.M.  Heimerl,  "Flat  Flame  Olympics:  Test  Problem  B,"  BRL  Technical  Report, 
to  be  published. 

Proceedings  to  be  published  in  Notes  on  Numerical  Fluid  Mechanics .  Vol .  5, 

Ed.  N •  Peters  and  J .  Wamatz ,  Vieweg-Verlag,  publication . 
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II.  THE  TEST  PROBLEM 


The  problem  is  an  unsteady  propagating  flame  with  one  step  chemistry 
(reactant  product).  The  Lewis  number  may  be  different  from  unity.  The 
governing  equations  are 


and 


3y  _  1  32Y 

3t  "  Le  .2 
dx 


9t 

9t 


+  R, 


(1) 


(2) 


where  Y  is  the  mass  fraction  of  the  reaction,  t  is  the  time  coordinate,  Le  is 
the  Lewis  number,  x  is  the  space  coordinate,  R  is  the  rate  of  the  reaction, 
and  T  is  the  temperature.  All  the  above  quantities  are  non-dimensional.  The 
reaction  rate  is  given  by 


R 


2Le 


Y  exp 


3  (1-T) 

1 -a ( 1 -T ) 


]  • 


The  boundary  conditions  (t>0)  are 


(3) 


T  =  0,  Y  =  1  at  x  =  -  00 


dT  dY 
dx  dx 


0  at  x  =  00  - 


(4) 


Since  we  are  interested  in  the  transient  problem,  initial  condition  must 
be  specified.  These  are  (t=0) 

Y  =  1  -  exp  (Le  x) 

T  =  exp  ( x )  ( 5 ) 

for  x  £  0  and 

Y  =  0,  t  =  1  (6) 

for  x  £  0. 

The  reaction  rate  was  normalized  such  that  the  steady  sj^ate  problem 
yields  the  flame  velocity  -1  in  the  asymptotic  limit  B  00 .  The  initial 


n 

W.B.  Bush  and  F.E.  Fendell,  "Asymptotic  Analysis  of  Laminar1  Flame  Propagation 
for  General  Lewis  Numbers,"  Combustion  Science  arid  TecfUj  Vol.  1, 
pp.  421-428,1970. 


8 


conditions  are  the  profiles  generated  by  the  asymptotic  analysis.  For  x  £0, 
it  is  assumed  that  the  chemistry  term  R  is  negligibly  small  (convective- 
diffusive  zone).  The  profiles  can  then  be  shown  to  be  exponential 
functions.  At  x  =  0,  it  is  assumed  that  adiabatic  conditions  have  been 
achieved.  There  are  discontinuities  in  the  first  derivatives  at  x  =  0. 

The  parameter  a  is  taken  to  be  0.8.  There  are  six  cases  to  be  solved; 
these  are  given  in  Table  1.  The  desired  output  is  the  flame  velocity/  the 
frequency  and  amplitude  of  flame  velocity  oscillations,  if  observed,  and  the 
enthalpy  profile.  In  the  above  normalization,  the  enthalpy  H  is  given  by 

H  =  Y  +  T.  (7) 

If  the  Lewis  number  is  equal  to  one,  we  can  add  Eqs.  (1)  and  (2)  and 
obtain 

9h  _  9jH  (8) 


Using  the  boundary  conditions,  we  can  obtain 

H  =  Y  +  T  =  1.  (9) 

So  the  Y  equation  can  be  eliminated. 


TABLE  1 .  THE  CHOICES  OF  THE  INPUT  PARAMETERS  TO  BE  USED 


Case 

B 

Lc 

1 

10 

1.0 

2 

10 

2.0 

3 

10 

0.5 

4 

20 

1.0 

5 

20 

2.0 

6 

20 

0.5 

Because  of  the  above  simplification,  the  choice  of  Lewis  number  equal  to 
one  is  a  common  assumption.  One  purpose  of  the  problem  as  given  is  to  observe 
what  effect  other  choices  of  the  Lewis  number  may  have. 


Q 


The  above  problem  is  somewhat  ambiguous  in  requesting  the  transient  flame 
velocity.  The  flame  velocity  is  defined  as  the  speed  at  which  the  flame  front 
propagates  with  respect  to  the  unburned  mixture.  At  steady  state,  this  is 
well  defined.  But  while  approaching  steady  state,  the  flame  front  can  change 
shape,  and  one  part  of  the  flame  travel  more  rapidly  than  another. 

We  arbitrarily  chose  to  compute  the  flame  speed  with  respect  to  T  — 

0.5.  This  is  the  average  of  the  unburned  and  the  adiabatic  temperature. 

It  will  be  useful  to  introduce  a  coordinate  transformation 


4'  -  x  +  x  (t) 
c 

xc  (0)  =  0. 


Then  by  the  chain  rule 

8  Y  8  Y  3i|j 
3t  +  3i|>  3t 

and 

8  t  3jr  8^ 

at  +  3i|>  3t 


i 


Le 


R 


+  R. 


(10) 


(ID 


(12) 


We  will  define  x  in  such  a  way  that  the  point  T  =  0.5  remains  at  the 
same  point  in  ^  space.  From  the  initial  conditions,  this  means  T  =  0.5  at 
=  in  (0.5).  The  details  of  this  procedure  are  given  in  the  next  section. 


We  will  use  the  notation 

8x 


SP(t)  = 


8\J/ 

3t 


8 1 


(13) 


The  equations  defining  the  flame  front  become 


3Y 
3 1 


SP 


8  Y 
8\|/ 


1 

Le 


82Y 


-  R 


and 

8t 
8 1 


+  R. 


(14) 


(15) 


SP(t)  is  the  speed  of  the  origin  of  the  (^/t)  coordinate  system  with 
respect  to  the  (x,t)  coordinate  system.  Since  the  flame  is  in  a  fixed 
location  in  the  (\|/,t)  system,  SP(t)  is  also  the  negative  of  the  velocity  of 
the  flame  (the  flame  propagates  toward  x  =  -  As  the  flame  approaches 

steady  state,  SP(t)  approaches  the  steady  state  speed  of  the  flame  front, 
and  ST/St  and  3y/St  approach  zero. 
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III.  THE  NUMERICAL  SCHEME 


The  package  PDECOL,  developed  by  Madsen  and  Sincovec,  was  used  to  solve 
the  above  equations.®  This  is  a  general  package  for  solving  partial 
differential  equations.  It  has  been  modified  so  as  to  efficiently  solve  flame 
problems . 


The  spatial  discretization  is  accomplished  by  finite  element  collocation 
methods  based  on  B-splines.^  The  basic  assumption  is  that  the  solution  can  be 
written  in  the  form 


Y 


NC 

I 


B,  (^)  and 
1 


T  = 


NC 

I 


(t)  B±  (T|;) 


(16) 


where  the  basis  functions  B± ( ^ ,  i  =  1,  ...NC,  are  B-splines  and  span  the 

solution  space  for  any  fixed  t  to  within  a  small  error  tolerance.  The  time 
dependent  coefficients  Cj^1^  are  determined  uniquely  by  requiring  that  the 
expansion  above  satisfy  the  given  boundary  conditions  and  that  they  satisfy 
Eqs.  (14)  and  (15)  exactly  at  (NC-2)  interior  (collocation)  points.  The  B^ 
are  piecewise  polynomials  of  order  KORD.  Given  a  user  supplied  set  of  NB 
breakpoints,  (i.e.,  a  set  of  strictly  increasing  locations  where  the 
polynomials  are  joined),  and  the  number  of  continuity  equations,  NCC,  to  be 
applied  at  the  breakpoints,  PDECOL  generates  a  set  of  NC  =  KORD  (NB-1)  - 
NCC(NB-2)  basis  functions  and  collocation  points.  Since  by  definition  a  B- 
spline  is  zero  except  over  a  small  interval,  at  any  collocation  point  no  more 
than  KORD  of  the t B-splines  are  non-zero.  So  the  system  of  ODE's  for  the 
coefficients  ^  will  not  be  fully  coupled. 


Fairly  general  boundary  conditions  are  allowed, 
into  ODE's. 


They  are  also  converted 


The  system  of  ODE's  is  integrated  in  time,  using  a  variant  of  the  Gear 
stiff  implicit  integrator.  The  appropriate  banded  Jacobian  is  generated 
internally  by  the  program. 


o 

N.K.  Madsen  and  R.F .  Sincovec ,  ” Algorithm  540 .  PDffCOL,  General  Collocation 
Software  for  Partial  Differential  Equations  [ D3]9  A  GMT  DM  Ss  Vol .  5,  pp .  326 - 
351,  1979 . 


^ 0 .  deBoor y  " Package  for  Calculating  with  B-Splines ;  Siam .  J .  Nurner*  Anal . 
Vol .  14,  pp.  441-472,  1977. 
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The  time  integration  is  controlled  by  a  user  supplied  error  tolerance 
Single  step  error  estimates  divided  by  CMAXk(^  will  be  kept  less  than 


e  in  the  root -mean-square  norm.  In  PDECOL,  CMAXk(l)  is  initially  set  to  the 
maximum  of  |ckU)|  and  1.0.  Thereafter,  CMAXk(l*  is  the  largest  value  of 
|ckl  ;|  seen  so  far,  or  the  initial  CMAXk(l)  if  that  is  larger.  This  error 
criterion  is  not  appropriate  for  flames.  In  some  locations  species 
concentrations  will  approach  zero,  and  we  do  not  want  to  waste  time  computing 
a  negligible  concentration  very  accurately. 


So  we  introduce  a  semi-relative  error  control.  CMAXk(i)  will  be  chosen 
as  the  maximum  of  ck(1'  and  a  user  supplied  parameter  SREC.  So  mass  fractions 
less  than  SREC  will  be  computed  less  accurately.  In  this  paper,  SREC  =  10~^ 

e  . 


The  original  program  PDECOL  is  fully  implicit.  That  is,  the  Jacobian  of 
the  set  of  ODE's  is  used  to  advance  the  time  integration.^ 

Our  main  interest  is  in  the  steady  state  solution  of  flames  with  detailed 
chemistry.  Such  flames  involve  relatively  large  numbers  of  chemical  species 
and  hence  PDE ' s .  As  the  number  of  PDE*s  increase,  both  the  overall  dimensions 
of  the  Jacobian  and  the  bandwidth  increase  rapidly.  This  creates  problems 
both  in  terms  of  storage  space  and  execution  time. 

So  we  essentially  uncouple  the  equations  and  solve  them  successively.  We 
assume  that  changes  in  one  PDE  do  not  effect  the  others.  That  is,  the  terms 
in  the  Jacobian  relating  to  the  effect  of  one  PDE  on  another  are  set  equal  to 
zero.  Then  the  Jacobian  matrix  can  be  decomposed  into  smaller  matrices,  one 
for  each  PDE. 

To  integrate  a  time  step,  the  predictor  of  the  predictor-corrector  method 
is  used  to  get  first  estimates  for  all  the  ck^1^  at  the  new  time.  Then  we 
solve  for  the  1  ,  using  the  first  small  Jacobian,  and  assuming  that  the 
other  ck  1  do  not  change.  Then  we  use  the  new  values  for  c.  the 

predicted  values  for  c^^,  k  >  2  ,  and  solve  for  c2^*t  This  continues  for 
all  the  PDE's.  The  error  is  estimated,  and,  if  necessary,  the  above  procedure 
is  iterated. 

This  procedure  is  similar  to  one  developed  by  Spalding  and  Stephenson  for 
use  in  a  finite  difference  code.^ 

For  the  present  problem,  there  are  only  two  PDE * s .  Moreover,  the  details 
of  the  transient  solution  are  required.  In  this  case  it  would  be  more 
efficient  to  use  the  full  Jacobian.  However,  we  did  not  want  to  rewrite  the 
code,  so  it  was  run  using  successive  calculation.  As  a  result,  smaller  time 
steps  probable  had  to  be  taken  to  obtain  the  same  accuracy. 


A.C.  Hindmarsh ,  " Preliminary  Documentation  of  GEARIB.  Solution  of  Implicit 
Systems  of  Ordinary  Differential  Equations  with  Banded  Jacobian,"  Rep.  UCID- 
30130,  Lawrence  Livermore  Laboratory,  1976. 


1:LD.B.  Spalding  and^P.L.  Stephenson,  "Laminar  Flame  Propagation  in  Hydrogen  + 
Bromide  Mixtures,  Proc.  R.  Soc.  Lond.  A .  Vol.  324,  pp.  315-337,  1971. 
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The  flame  front  will  be  kept  near  the  center  of  the  interval  of 
integration,  so  the  breakpoints  should  be  densest  here.  Rather  then  choosinq 
the  breakpoint  sequence  by  hand,  an  algorithm  has  been  written  to  generate  an 
appropriate  set  of  breakpoints. 

The  user  supplies  L,  the  total  length  of  the  interval  of  integration; 
NINT,  the  total  number  of  subintervals  (NB=NINT+1);  NCN,  the  number  of 
subintervals  of  equal  length  that  will  be  at  the  center  of  the  interval  of 
integration;  and  FC,  the  ratio  between  the  longest  subintervals  (on  the 
boundaries)  and  the  shortest  subintervals. 

The  program  generates  a  set  of  subintervals  whose  lengths  increase  by  a 
constant  factor  Y ,  where 

Y  =  log-1  ( 2 ( log  FC)/ (NINT-NCN) ) •  (17) 

The  commmon  length  LC  of  the  NCN  shortest  subintervals  in  the  center  is 

.  .  ^  ,  ( NINT-NCN )  /2 

LC  =  T  /  (NCN  +  2Y(Y  -1)/(Y“D).  (18) 

To  actually  run  the  code  for  the  transient  problem,  the  user  chooses  the 
defining  parameters  a, 3,  and  Le.  Then  the  numerical  parameters  L,  NINT,  NCN, 
FC,  KORD,  NCC ,  and  TFINAL  must  be  chosen.  TFINAL  is  the  time  t  the  code  will 
integrate  to.  In  this  paper  we  always  choose  K0RD=4  and  NCC  =  2  (cubic 
splines) . 

The  boundary  values  and  are  chosen  as  approximately  -L/2  and 

L/2.  These  are  adjusted  slightly  so  the  value  £n  (.5)  is  a  collocation 
point.  So  the  value  T  =  0.5  is  at  a  point  where  the  computation  is  most 
accurate • 

As  the  time  integration  proceeds,  we  require  that  3T/3t  equal  zero  at 
this  collocation  point.  This  is  done  by  adjusting  the  value  of  SP.  At  t  =  0, 
we  choose  SP(t)  =  0,  then  integrate  one  time  step.  At  the  end  of  the  time 
step,  evaluate  3T/3t  and  3t/3^,  using  the  expansion  (16)  and  the  space 
derivative  of  the  expansion. 

Now  choose 

SPF  =  ( 3 t/3 t )  /  (3T/3^).  (19) 

If  SP  was  given  this  value  at  t  =  0,  the  present  time  derivative  would  be 
zero.  SPF  is  taken  as  the  flame  speed  at  the  present  time.  At  later  time 
steps  increments  to  SPF  are  computed  in  the  same  fashion. 

The  flame  speed  is  then  computed  with  a  lag  of  one  time  step.  It  is 
possible  for  the  lag  to  cause  the  value  T  =  0.5  to  drift  away  from  a  fixed 
location. 

So  a  heuristic  correction  factor  is  used.  We  define 

FVD  =  1  0 ( . 5  -  T) /TFINAL  (20) 


and  let 
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SP  =  ( 3 t/3 t  -  FVD)/(9t/9^) 


(21) 


This  value  is  actually  used  in  the  equations. 

_*3 

Even  with  a  moderate  error  control  (£  =  10  )  SPF  and  SP  agree  to  3  or  4 

places.  So  the  drift  of  the  value  T  =  0.5  due  to  the  lag  in  computing  SP  has 
a  negligible  effect. 


This  procedure  can  be  checked  at  steady  state 
Integrating  Eq.  (14),  one  obtains 


SP  (Y| 


-  Y 


R 


—  i 

94) 


R 


-  —  I 

94<  '4) 


( 9  Y/9 1  =  9T/9t  =  0). 

4> 

Jt  R  Rd4>,  (22) 


The  integral  is  evaluated  using  the  trapezoid  rule.  At  steady  state,  the 
value  obtained  by  this  procedure  should  be  the  same  as  the  value  computed  from 
the  T  equation. 


IV.  NUMERICAL  SOLUTIONS 

The  test  problems  were  solved  using  a  Cyber  76.  The  input  parameters 
used  are  given  in  Table  2.  Table  3  gives  the  run  times,  the  number  of  time 
steps  taken,  and  the  steady  state  flame  speeds.  The  speeds  are  close  to  the 
asymptotic  limit  of  1 .  As  expected,  they  are  closer  to  1  for  B  =  20  then  for 
B  =  10. 

Numerical  errors  can  be  introduced  because  of  inadequate  spatial 
resolution  (NINT,  NCN,  and  FC)  or  inadequate  temporal  resolution  (£).  Errors 
can  also  occur  if  L  is  assigned  a  value  that  is  too  small.  The  profiles  may 
then  have  nontrivial  space  derivatives  at  the  boundaries.  The  boundary 
conditions  will  then  distort  the  flame  front. 

As  a  check  on  the  numerical  accuracy,  the  problems  were  rerun  with  NINT 
decreased  by  25%,  £  increased  by  a  factor  of  3,  and  L  decreased  by  10%.  The 
results  were  virtually  identical  to  those  reported  here. 

Figure  1  shows  the  flame  speed  profile  for  Case  1.  The  profile  is  only 
shown  up  to  t  =  0.5.  After  this  time  there  is  just  a  slow  approach  to  the 
final  steady-state  value. 
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THE  COMPUTATIONAL  INPUT  PARAMETERS  FOR  THE  TEST  PROBLEMS 


NINT 

e 

L 

NON 

FC 

TF INAL 

40 

i 

o 

20 

20 

10 

10 

40 

1 

O 

20 

20 

10 

50 

60 

10-4 

40 

30 

30 

100 

40 

1 

O 

T— 

20 

20 

20 

10 

80 

1 

o 

24 

40 

40 

100 

60 

1 

o 

r— 

40 

30 

40 

50 

TABLE 

3.  RESULTS 

FOR  THE 

TEST  PROBLEMS 

Run  Time 

No.  of 

Steps 

SP 

24.1 

568 

.918 

18.5 

500 

.849 

60.0 

898 

.958 

18.7 

478 

.954 

404.5 

6597 

63.7 

898 

.978 

IS 


Figure  2  shows  the  temperature  profile  for  Case  1  for  selected  values  of 
t.  At  t  =  0  there  is  a  slight  Gibbs  overshoot  at  x  =  0  where  the  exponential 
function  is  joined  to  the  straight  line.  The  steep  gradient  near  x  =  0  causes 
an  increase  in  the  flame  speed  due  to  thermal  conductivity.  This  steep 
gradient  relaxes  rapidly.  There,  there  is  a  relatively  slow  approach  to  the 
steady  state  solution.  This  is  quite  close  to  the  asymptotic  solution  except 
near  x  =  0. 

For  Case  1  (Le  =1  ),  Y  =  1  -  T  and  H  =  1,  so  these  profiles  are  not 
shown. 

For  Case  2  (Le  =  2),  after  the  initial  flame  speed  increase  the  speed 
drops  off  rapidly  (Figure  3).  It  then  takes  a  relatively  long  time  to 
approach  steady  state  (Figure  4) . 

This  occurs  because  the  diffusion  of  the  reactant  is  now  slow.  After  the 
initial  burst  in  the  flame  speed  the  reactant  becomes  used  up  near  x  =  0  and 
the  flame  front  has  not  yet  moved  very  far.  So  the  temperature  in  the  flame 
front  decreases.  Finally,  the  flame  propagates  forward,  more  Y  becomes 
available,  and  the  temperature  profile  increases  to  steady  state  (Figures  5 
and  6).  The  enthalpy  profile  (Figure  7)  behaves  similarly.  This  kind  of 
overshoot  can  only  occur  because  the  Lewis  number  is  different  from  unity. 

For  Case  3,  the  flame  speed  profile  (figure  8)  is  similar  to  Case  1. 
Diffusion  is  now  faster  than  thermal  conductivity,  and  there  is  no  overshoot 
due  to  lack  of  reactant. 

It  is  necessary  to  increase  L,  the  length  of  the  interval  of 
integration.  Since  diffusion  is  more  rapid,  the  Y  profile  is  more  spread 
out.  Figure  9  shows  the  Y  profile  over  the  same  interval  as  for  Case  2. 
However,  the  numerical  interval  must  extend  further  to  resolve  the  approach  to 
the  cold  boundary  condition  Y  =  1.  Here  we  also  increase  FC,  the  ratio 
between  the  center  and  boundary  subintervals,  so  that  the  flame  front  is  still 
adequately  resolved. 

Figure  10  shows  the  T  profile,  which  is  not  as  spread  out.  The  enthalpy 
profile  is  given  by  Figure  11. 

The  next  three  cases  (3  =20)  have  a  higher  activation  energy,  so  the 

asymptotic  analysis  should  be  more  accurate.  This  is  true  for  Case  4  (Le  = 

1).  The  flame  speed  profile  shows  only  a  small  excursion  before  reaching 
steady  state  (Figure  12).  The  final  T  profile  (Figure  13)  is  almost  identical 
to  the  asymptotic  solution. 

For  case  5  (Le  =  2)  the  flame  speed  profile  begins  similarly  to  Case  2 
(Figure  14).  However,  the  flame  velocity  then  oscillates  (Figure  15)  with  a 
period  of  about  9  time  units . 

This  effect  is  due  to  the  slowness  of  the  diffusion  of  Y  combined  with 
the  high  activation  energy  of  R.  The  Y  and  T  profiles  first  spread  out  due  to 
diffusion  and  thermal  conductivity.  Then  the  reaction  cuts  in.  Due  to  the 
rapidity  of  the  reaction,  a  bump  occurs  in  the  T  profile.  In  this  high 
temperature  region,  the  supply  of  Y  is  exhausted.  So  the  temperature  profile 
again  relaxes  (Figures  16  and  17)  and  the  profiles  oscillate  around  the 
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asymptotic  solution  but  do  not  converge  to  it.  The  enthalpy  profile  (figure 
18)  behaves  similarly. 

The  last  case  again  converges.  The  Y  profile  is  spread  out  by  diffusion, 
but  the  numerical  solution  is  again  very  similar  to  the  asymptotic  solution. 

V.  CONCLUSIONS 

The  results  show  that  for  one  step  chemistry  and  high  activation  energy 
the  asymptotic  analysis  is  very  accurate.  Since  the  asymptotic  analysis  does 
not  resolve  the  flame  front,  the  steep  gradient  does  cause  a  transient 
increase  in  the  flame  speed  before  the  steady  state  solution  is  achieved. 

If  Le  <  1,  no  additional  problems  occur,  except  that  convergence  is 
somewhat  slower.  However,  the  combination  of  Le  >  1  and  a  high  activation 
energy  can  cause  overshoot  or  even  large  oscillations. 
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Figure  1.  Case  1.  Flame  Speed  profile. 
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Figure  2.  Case  1.  Temperature  Profile.  Time  =  0  (line).  Time 
0.032  (dot).  Time  =  10  (chain-dot). 
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Figure  3.  Case  2.  Flame  Speed  Profile  (0  <  time 
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Figure  5.  Case  2.  Reactant  Profile.  Time  =  0  (line).  Time 
0.029  (dot).  Time  =  0.995  (chain-dot).  Time  = 

50  (dash). 
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Figure  6.  Case  2.  Temperature  Profile.  Time  =  0 

(line).  Time  =  0.029  (dot).  Time  =  0.995  (chain- 
dot)  •  Time  =  50  (dash) . 
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Figure  7.  Case  2.  Enthalpy  Profile.  Time  = 

0  (line).  Time  =  0.029  (dot).  Time  =  0.995  (chain- 
dot).  Time  =  50  (dash). 
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Figure  8.  Case  3.  Flame  Speed  Profile. 
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Figure  9,  Case  3.  Reactant  Profile.  Time  =  0  (line).  Time 
0.042  (dot).  Time  =  50  (chain-dot). 
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Figure  10.  Case  3.  Temperature  Profile.  Time  =  0  (line).  Time 
=  0.042  (dot).  Time  =  50  (chain-dot). 
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Figure  11.  Case  3.  Enthalpy  Profile.  Time  =  0  (line).  Time 
0.042  (dot).  Time  =  50  (chain-dot). 
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Figure  12.  Case  4.  Flame  Speed  Profile. 
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Figure  13.  Case  4.  Temperature  Profile.  Time  =  0  (line) .  Time 
0.034  (dot).  Time  =  10  (chain-dot). 
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Figure  14.  Case  5.  Flame  Speed  Profile  (0  <  time 
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Figure  15.  Case  5.  Flame  Speed  Profile  (0  <  time  <  100). 
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Figure  16.  Case  5.  Reactant  Profile.  Time  =  0  (line).  Time 
4.28  (dot).  Time  =  8.98  (chain-dot).  Time  = 

13.483  (dash). 
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Figure  17.  Case  5.  Temperature  Profile.  Time  =  0  (line).  Time 
=  4.28  (dot).  Time  =  8.98  (chain-dot).  Time  =  13.483 
( dash) . 
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Figure  18.  Case  5.  Enthalpy  Profile.  Time  =  0  (line).  Time  = 
4.28  (dot).  Time  =  8.98  (chain-dot).  Time  =  13.483 
( dash) . 
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Figure  19*  Case  6.  Flame  Speed  Profile 
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Figure  20.  Case  6.  Reactant  Profile.  Time  -  0  (line).  Time 
0.045  (dot).  Time  =  50  (chain-dot). 
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Figure  21.  Case  6.  Temperature  Profile.  Time  =  0  (line).  Time 
=  0.045  (dot).  Time  =  50  (chain-dot). 
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Figure  22.  Case  6.  Enthalpy  Profile.  Time  =  0  (line).  Time 
0.045  (dot).  T.^~~  =  50  (chain-dot). 
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